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Vortices carrying fractions of a flux quantum are predicted to exist in multiband superconductors, 
where vortex core can split between multiple band-specific components of the superconducting con¬ 
densate. Using the two-component Ginzburg-Landau model, we examine such vortex configurations 
in a two-band superconducting slab in parallel magnetic field. The fractional vortices appear due 
to the band-selective vortex penetration caused by different thresholds for vortex entry within each 
band-condensate, and stabilize near the edges of the sample. We show that the resulting fractional 
vortex configurations leave distinct fingerprints in the static measurements of the magnetization, as 
well as in ac dynamic measurements of the magnetic susceptibility, both of which can be readily used 
for the detection of these fascinating vortex states in several existing multiband superconductors. 

PACS numbers: 74.25Ha, 74.25Uv, 74.70Xa, 74.70Ad 


Multiband superconductors 0,0 present a variety of 
intriguing properties that are not found in their single¬ 
component counterparts. Theoretical predictions have 
added more striking properties to that list and challenge 
experiments to prove them. One of such properties is the 
appearance of fractional vortices in multiband materials 
|3(, seemingly violating flux quantization. This is only 
possible for different winding numbers of different order 
parameters in a system of coexisting weakly interacting 
condensates, and is facilitated for significantly different 
length scales of the condensates - especially under meso¬ 
scopic confinement 0 - 0 - Weakly coupled multiband ma¬ 
terials 0 and superconducting multilayers [llj as their 
artificial analogue are readily available, hence clever ex¬ 
periments should be devised for detec ting and manipu¬ 
lating fractional vortices (see e.g. Ref. 1121 ). In addition, 
dynamic dissociation of vortices is predicted in the flux 
flow regime [13j |, as well as the stationary vortex splitting 
00 stemming from phase frustration in superconduc¬ 
tors with three or more bands 0 0, but neither of 
those vortex fractionalizations has been realized to date. 


In this Letter we explore the effect of a surface in stabi- 
lizin g th e fractional vortices in multiband superconduc¬ 


tors [lSj, and propose static ( dc ) and dynamic (ac) mea¬ 
surements to directly detect them. We consider a two- 
band superconducting slab in parallel magnetic field H , 
with width much larger than the field penetration depth 
in order to prevent strong confinement effects. For our 
numerical experiments, we have used the two-component 
Ginzburg-Landau (TCGL) model, where by cautiously 
setting temperature T close to the critical temperature 
T c , we ensure the qualitative and quantitative validity of 
our results (for comparison with other available theoret¬ 


ical models, see e.g. Re fs. 0 - 0 ). In the TCGL frame¬ 
work, as given in Ref. l22l eight independent material 
parameters are needed for a system with both interband 
and magnetic coupling, namely, the Fermi velocity of the 
first band iq, the square of the ratio of the Fermi ve¬ 
locities in the two bands a = ( — ) 2 , the elements of the 
coupling matrix An, A22 and A12 = A21, the total density 
of states N{ 0) as well as the partial density of states of 
the first band n\ (712 = 1 — n 1 ), and finally T c , which sets 
the energy scale W 2 = 8n 2 T//7((3). The TCGL free 
energy functional reads 0] 

i= 1,2 3 

—r(^rv , 2 + ^1^2) + ( h i (i) 

where j = 1,2 is the band index, ay = —N(0)rijXj = 
—N(0)rij(T — Sj/rij6), /3j = ( N(0)rij)/W 2 , rrij = 
3W 2 /(N(0)rijV 2 ) , and V = (iV(0)Ai2)/<5, with S being 
the determinant of the coupling matrix, and S, Si and 
S 2 defined as in Ref. 0. The local magnetic field in the 
sample is denoted by h and the external applied field by 
H. 

Minimization of the free energy in Eq. © with re¬ 
spect to t/jj and A yields the Ginzburg-Landau equations. 
Introducing the normalization for the order parameters 
by W, for the vector potential by A 0 = hc/Aeir/i, for 
the lengths by £1 = hvi/^GW, and for the time by 
to = 4 tt(jkI(i/c 2 (cr is the normal-state conductivity), 
the dimensionless time-dependent TCGL equations in 
the zero-electrostatic potential gauge are written as: 
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V^- = (-tV - A)Vi - (XI - IV’il 2 )^! - 7^2, (2) 

= ~{-iS7-A) 2 *l) 2-(x2-|'02| 2 )'i/’2--^y%^i, (3) 

at a ar/sj 

dA 

-fo = 3°~ K i V x V x A, (4) 

where = 7^^ 2ni N(o) ’ K 2 = K 1 a y /n 1 /n 2 , 7 = 

Ai2/ni5, and 77 = Trh/(8t 0 T c ). In Eq. (j4j) the super¬ 
current density is, 



; s = K[il> 1 (z V - 1)^1*] + ^^ 2 (z V - A)r 2 \ , (5) 


where 7?. denotes the real part of the expression. After 
the made choice of normalization units, we are left with 
seven parameters: An, A 22 , A 12 , vi/v 2 , ni, N( 0) and 
r). We fixed T = 0.85T C to firmly remain in the validity 
regime of the TCGL theory. For the other parameters, 
we take An = 2.0, A 22 = 1.03, A 12 = 0.005, v\/v 2 = 0.52, 
rii = 0.355 and 77 = 5.0, while N( 0) is fixed by chosen 

Kl = 10.0. 

In our numerical experiment, we studied a supercon¬ 
ducting slab of width 100Ci, corresponding to 7.42A for 
the considered parameters (A is the magnetic penetration 
depth), in the presence of a parallel time-dependent mag¬ 
netic field H(t) = Hdc + Hac cos(uit) (with frequency unit 
wo = 1/io)- The TCGL equations j2])-(|4|) were integrated 
on a two dimensional grid with grid spacing a x = a y = C, 1 , 
much smaller than any characteristic length scale at the 
considered temperature. The discretization was imple¬ 
mented by the link variable method which preserves the 
gauge invariance of these equations Q- For the iter¬ 
ative solver, we combined a relaxation method with a 
stable and accurate semi-implicit algorithm J25|. Peri¬ 
odic boundary conditions were applied in the x direction 
(with size of the unit cell 200£i) whereas for the y direc¬ 
tion we imposed Neumann boundary conditions at the 
superconductor-vacuum interface (for details of the nu¬ 
merical implementation, please see Ref. E3) . The sub¬ 
sequently calculated magnetization, M = ((h) — H)/4n 
((...) denotes spatial averaging inside the sample), is a 
measure of the expelled flux from the sample and the cor¬ 
responding M(H) response was obtained by ramping up 
the magnetic field with steps of AH = 2 x 10 ~ 4 (in units 
of Hq = hc/2eC, 2 ). For the study of magnetic relaxation 
dynamics (for H ac ^ 0), we calculated the imaginary 
part of the magnetic susceptibility as the Fourier trans¬ 
form of M(t), x"(Hdc,u ) = -^A- / 0 27r M(t)sm(ut)d(ut). 
x" is directly proportional to the time average of the 
energy dissipated in the sample, as can be seen from 
the expression for the energy dissipated in one cycle, 
W = 4n f mcLH oc x" ■ The local dissipation of energy, 
W(R,t), comprises two terms 26], 271, one is the Joule 


FIG. 1. (Color online) Magnetization versus applied magnetic 
field for a two-band superconducting slab, at T = 0.85T C , 
comprising stability curves of the obtained different vortex 
configurations (dotted line shows the sequence of states in 
increasing magnetic field). The red curves correspond to the 
fractional vortex states whereas the blue ones correspond to 
the composite vortex states. Labels (a)-(t) are used to denote 
different vortex states. 


heating term due to the normal currents, proportional to 
\dA(R, t)/dt\ 2 , and the other is related to the relaxation 
of the order parameter, proportional to \dip(R,t)/dt | 2 . 

The calculated M(Hd c ) (for H ac = 0) in units of the 
thermodynamic critical field H c is shown in Fig. [T] 
which for increasing magnetic field follows the dotted 
line. It exhibits a series of steps corresponding to the 
entry of fractional vortices, and forming vortex configu¬ 
rations shown in Fig. [2j The field for the first vortex 
penetration H p = 0.764 H c is superheated due to surface 
effects, and, as the coherence lengths associated with the 
two band-condensates differ from each other significantly 
(£2 = 2.24£i) for the here considered parameters, the 
vortex entry first occurs in the second band-condensate, 
where surface barrier is suppressed at a lower magnetic 
field. Consequently, the vortex configuration after the 
first jump in M ( Hd c ) consists only of fractional vortices 
in the second band-condensate (as shown in Fig. [2])a2) 
). The fractional vortices find their equilibrium positions 
near the surface in a similar fashion to those reported 
in Ref. Il8l . where the London theory was used in the 
absence of interband coupling. However, the here cal¬ 
culated penetration field H p (= 0.764 H c ) is larger than 
the one predicted in Ref. 18 (0.521 H c ), which is expected 
since London approach neglects the influence of the finite 
size of the vortex cores and hence does not capture the 
corresponding energy needed for the vortex entry. As the 
magnetic field is further increased, the fractional vortices 
of the first band-condensate also penetrate the sample 
(at Hdc = 0.772.ff c ), and combine with those of the sec¬ 
ond band to form composite vortices, which afterwards 
further penetrate the central part of the sample (see Fig. 
0(b) ). Similar scenario continues at higher magnetic 
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FIG. 2. (Color online) Vortex configurations corresponding 
to selected states from Fig. [T] For each state, left/right panel 
shows the Cooper-pair density of the first/second band con¬ 
densate, and are respectively tagged 1 and 2. 


field, where fractional and composite vortex states are 
alternately stabilized in the sample, as shown in Figs. |T| 
and[2j Beyond state p in Fig. |I] the second-band con¬ 
densate is fully depleted (see Fig. [2](q2) ), and fractional 
vortex states are no longer possible. 

In fact, one can directly obtain quantitative informa¬ 
tion about the fractional flux states from the magneti¬ 
zation curves shown in Fig. [TJ To do so, we first es¬ 
timated (for the considered parameters) the fraction of 
the flux quantum <f>o carried by the vortex in each band- 



FIG. 3. (Color online) (a) Graph illustrating the calculated 
difference —AtyAM, taken between the magnetization curve 
of state a and the reconstructed Meissner curve of Fig. Q] (b) 
Graph showing direct link between the calculated —AttAM 
(related to the flux entry into the sample) and the distance 
d of the fractional vortices to the surface (extracted from the 
respective vortex configurations), along the stability curve of 
state a. 


condensate, i.e. = 0.28<f>o and $2 = 0.724>o, and 
compared those fractions with the total flux entering the 
sample at the point of nucleation of e.g., the state of 
Fig. [2](a), estimated from the jump in magnetization. 
Note that the Meissner curve shown in Fig. [T] exhibits a 
nonlinear behavior near the penetration field due to the 
depreciation of the order parameter near the edges [28j . 
so, to obtain the correct estimate of flux entry one needs 
to consider M oc Hd c for the entire Meissner curve. The 
entering magnetic flux is then calculated from the differ¬ 
ence in magnetization AM between the (reconstructed) 
Meissner curve and the a branch of Fig. |T] at the flux 
penetration field, as shown in Fig. [3ja) , amounting to 
4.4<f>o- Since Fig. [2](a2) shows 8 penetrating fractional 
vortices in the second band-condensate, we obtain mag¬ 
netic flux per fractional vortex of (f >2 = 0.55$o - which 
is lower than our first estimate of $2 = 0.72. Such a 
reduced value for </>2 is related to the proximity of the 
fractional vortex to the surface, and the interaction of 
its current with the screening currents running at sur¬ 
faces. Introducing the correction due to screening cur¬ 
rents <f >2 = 4> 2 [1 — exp(—d/A)], with d the vortex distance 
from the surface extracted from the calculated vortex 
configuration (see Ref. [29h . we obtain $2 « 0.7$o in the 
entire magnetic field range of stability for state a , which 
compares very well with our first estimation of $ 2 - The 
vortex proximity to the surface also explains why AM 
shown in Fig. G^a) diminishes for decreasing Hd c , as the 
fractional vortices approach the surface for lowered field 
and the flux carried by them also decreases. As we show 
in Fig. [31(b), distance d is approximately linear with Hdc , 
and proportional to AM. Therefore, based on this un¬ 
derstanding even the location of the fractional vortices 
in the sample can be deduced from the static magneti- 
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FIG. 4. (Color online) (a) Magnetic susceptibility x" as a 
function of the frequency u> of the ac magnetic field, calculated 
for composite vortex state and the fractional vortex state (for 
four equidistantly taken values of Hdc along the b (0.797 H c < 
Hdc < 0.821 H c ) and c (0.837 H c < Hdc < 0.886 H c ) branches 
in Fig. [l] in both sets of curves, x" was larger for larger Hdc)- 
(b) x" for fixed lo (several indicated values), as a function of 
the dc magnetic field. In both panels, the amplitude of the 
ac magnetic field was kept constant at H ac = 0.008 H c . 


zation data, provided that the number of participating 
fractional vortices is known from the start. For ideal 
surfaces, that number will correspond to the number of 
flux quanta deduced from the difference in magnetization 
between the reconstructed Meissner curve and the first 
composite vortex state (in the present case, state 6, see 
Fig. [Da) ). 

In what follows, we turn to the study of the ac mag¬ 
netic response of the vortex configurations found in the 
M(Hdc) curve, via calculations of x"(w, Hd r ). Our dis¬ 
cussion here is based on Refs. (33 and (33 that studied 
theoretically the ac dissipation in single-band mesoscopic 
superconductors. For such a system, the main contribu¬ 
tion to ac losses comes from the vortex nucleation re¬ 
gions near the sample edges and not from the vortices 
located inside the sample. As shown in Ref. 133, the 
frequency dependence of presents two peaks for a 
fixed Hdc near the vortex penetration field. One peak 
is found near ui = uio = c 2 /AnaX 2 and is related to the 
normal/superconducting current oscillations and conse¬ 
quent Joule heating. This peak is well described by the 
two-fluid model of superconductivity, and appears at fre¬ 
quency given by the characteristic time for the relaxation 
of the magnetic vector potential (see analysis in Ref. -30). 
Another peak, located at frequency oj p < wo, is due to 
the irreversible variation of the condensate wave func¬ 
tion during ac oscillations, since dissipation is intimately 


connected to the intrinsic relaxation time of the super¬ 
conducting condensate towards equilibrium [31|. Conse¬ 
quently, the frequency uj p depends on the relaxation rate 
of the order parameter (rj in Eqs. ([2E]), which is directly 
proportional to the Fermi velocity squared, and inversely 
to normal-state conductivity, a , and critical temperature 
T c ). It is also worth mentioning that the spatial distri¬ 
bution of losses inside the single-band mesoscopic super¬ 
conductor, W(R, t), shows that the lower the local order 
parameter near the surface, the larger is the dissipation 
corresponding to that region. 

Bearing above in mind, we performed similar analysis 
for our two-band slab. We take the ac magnetic field H ac 
much smaller than Hdc , ensuring that the system is in the 
linear regime of magnetic excitation. In our case, we still 
expect that the main contribution in x" arises from the 
vortex nucleation area near the surfaces, as argued above 
for the case of mesoscopic samples, since our sample has 
a moderate surface-to-volume ratio. On the other hand, 
our sample is large enough to allow entering fractional 
vortices to remain close to the edges, hence strongly in¬ 
teracting with the dissipative area where superconduct¬ 
ing condensate is depleted. Due to latter, we anticipated 
clear signatures of fractional vortices in the \"(u;,FOc) 
response. To prove this, we chose two branches in the 
M(Hdc) loops: one corresponding to the composite vor¬ 
tex state b and another to fractional vortex state c in Fig. 
[I] For fixed H ac = 0.008 H c and Hd c in the stability range 
of the studied states, we varied w and calculated x" ■ As 
shown in Fig. ID) a), by changing the frequency of the 
ac magnetic field, w, the composite vortex state always 
presents two dissipation peaks whereas for the fractional 
state the first peak is washed out. The second peak is 
common to all states since it originates from the nor¬ 
mal/superconducting current oscillations in which both 
band-condensates contribute, it does not depend on the 
number of vortices in the system, and it is even present 
in the Meissner state. 


The peak at lower frequencies is pronounced for the 
composite states due to the large area of suppressed 
second-band condensate (\ip 21 2 ) at the sample edges, sig¬ 
nalling large dissipation for further penetrating vortices. 
Actually, for the parameters we took, the dissipative be¬ 
havior for both composite and fractional vortex states 
will be governed by the second-band order parameter be¬ 
cause of its large coherence length, since the second-band 
condensate is more susceptible to the magnetic field than 
the first-band condensate. Thus, the difference in magni¬ 
tude of x"(w) (near its first peak) between these two con¬ 
sidered vortex states is due to the difference in depletion 
of the second-band condensate at the surfaces. For the 
composite states, where the two-peak structure in x"(ca) 
is pronounced, the screening supercurrents strongly de¬ 
plete the second-band order parameter at the surfaces, 
causing large dissipation near the first peak of x"(w). On 
the other hand, when the fractional vortices are stabilized 
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close to surfaces, the screening supercurrent diminishes 
there, the second-band order parameter is less depleted 
(c.f. Fig. (2Ta2) and Fig. [Hb2) ), and consequently, the 
dissipation peak is reduced. 

The difference in the dissipation of the two kinds of 
vortex states is even more evident for fixed u and varied 
Hdci which is a more suited experimental procedure. In 
this case, shown in Fig. [4jb) , we observe a sequence 
of peaks and valleys, following the exact sequence of 
composite and fractional vortex states from the M ( Hdc ) 
curve in Fig. [T] For composite states, x"(Hdc) rises with 
magnetic field, indicating high dissipation due to the in¬ 
creasingly depleted second-band condensate near the sur¬ 
face. Upon the penetration of the second-band vortices, 
fractional state is formed, x"(Hd c ) abruptly drops, and 
shows weaker dependence on Hd c - As a result, a remark¬ 
able profile of alternating peaks and valleys is obtained in 
x"(Hdc), very different from the simpler saw-tooth profile 
characteristic of mesoscopic single-band superconductors 

mu. 

In summary, we calculated static and dynamic mag¬ 
netic response of a two-band superconducting slab, and 
reported distinct properties that can be used to detect 
the fractional vortex states in multiband superconduc¬ 
tors. In static magnetometry, we showed how the analy¬ 
sis of the observed jumps in magnetization can be used to 
determine the fractional flux carried by vortices in differ¬ 
ent band-condensates, and the location of the fractional 
vortices with respect to the sample edge. Introducing an 
ac perturbation to external magnetic field, we demon¬ 
strated that the imaginary part of the magnetic suscep¬ 
tibility can identify fractional vortex states, both in its 
dependence on ac frequency and on dc magnetic field. 
Considering that recent superconducting materials are 
predominantly multiband (metal borides, iron pnictides, 
chalcogenides, etc.), our findings will stimulate further 
efforts in detection, manipulation and understanding of 
vortex states, creep and dynamics in those materials, as 
a precursor to potential applications. 
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